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appears in stochastic processes. 



cond-mat/ 9706 114 
June 1997 



* work supported by the DAAD programme HSP II-AUFE 

work done with partial support of the EC TMR programme, grant FMRX-CT96-0012 



1 



1. Introduction 

Several years ago Krug JTJ suggested the existence of boundary induced phase transitions 
in one-dimensional steady states. A line of first-order phase transitions (the so-called 
coexistence line) was found in the one-species asymmetric exclusion model P, [| |J] in 
which particles (call them positive) hop among vacancies (call them negative particles). 
On the coexistence line, the system is CP-symmetric (C corresponds to changing the 
sign of the particles, P is the parity operation) and the symmetry is spontaneously 
broken (some aspects of this phase transition which are relevant to this paper are 
reviewed in Appendix A). Since we want to compare the nature of the phase transitions 
in equilibrium and non-equilibrium phenomena, the phase transition just described 
would correspond to the two-dimensional n-state Potts model || with n > 4 at the 
critical temperature. The exclusion model does not have the equivalent of the low 
temperature domain which corresponds to the broken phase. 

The two-species model of Evans, Foster, Godreche and Mukamel || [7] presents 
a broken phase (the equivalent of the low temperature domain of the n-states Potts 
models) and this motivated us to study this model in detail. In our investigation we got 
a few surprises. 

We will just reproduce from Ref. |7|] the definition of the model, introducing also 
symmetry breaking boundary terms, and send the reader to the same paper in order 
to find out why the model is physically relevant. Each site of a one-dimensional lattice 
with L sites may be occupied by a positive particle or a negative particle or be empty. In 
each infinitesimal time step dt the following events may occur at each nearest-neighbor 
sites k, k + 1: 



all with probability dt. Where (+)&, (— )k an d (0)^ indicate a positive particle, a 
negative particle or a vacancy at the site k. In the same time step dt the following 
events may occur at the left boundary (k = 1) 



(+)*(0) fc+ i 



(0) fc (+)fc+i 




(1.1) 





(+)i with probability a dt 




(1.2) 



and at the right boundary (k = L) 



( )l with probability a dt 

->• ( ) L with probability (3(1 + h) dt . 



(1.3) 



One notices that for h — the probabilities are CP-invariant, the case h bigger than 
zero will be considered in Section 4. 

The time evolution of the system is given by a master equation or its equivalent 
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imaginary time Schrodinger equation 
d 



dt 



P >= -H\P > 



;i-4) 



where \P > is the probability vector and H is the hamiltonian of a one-dimensional 
quantum chain which is given in Appendix B. For most of this paper we will be interested 
in the properties of the steady state. 

Let us introduce some notations: p(k), m(k) and v(k) will denote the steady state 
density of positive, negative particles and vacancies at the site k. Their average values 
being p, m respectively v. It is also useful to define the quantities 

d(k) = p(k) — m(k) 

s(k) = (p(k)+m(k))/2, (1.5) 
their average values being d and s. 




Figure 1. Phase diagram for h=0. The phases A, B, C and D are obtained in 
mean- held and dehned in the text. Since the C phase is very narrow it is shown in an 
insert. The diamonds describe the boundary between the phase B and the phase D as 
obtained in Monte Carlo simulations. (No phase C is observed.) 

The phase diagram corresponding to the steady state was obtained in mean-field 
theory for the CP-symmetric case and is shown in Fig.|l|. One distinguishes four 
phases, we will give the average densities of particles in each one. The power law phase 
(A): 

p = m = -, v = (1.6) 
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The low densities (p = m < 1/2) symmetric phase (B): 

V = m = ——7 { <~ (1-7) 
a + p 2 

The low densities (p, m < 1/2) broken phase (C): Here mean-field gives three solutions, 
A symmetric unstable solution 

p = m = a ^ — < - (unstable), (1.8) 
a + p 2 

and two stable solutions 

a[3 



2 1- 



d = ±2 




1.9) 



The low density/high density phase (D): mean- field gives again three solutions. A 
symmetric unstable one which is again described by equation ( |1.8[ ) like in phase C, a 
second solution in which the positive particles have a high density (p > 1/2) and the 
negative particles have a low density (m < 1/2): 



p=l-P 

1 + a 1 

m = 

2 2 

and a third solution in which in Eq. ( |1.10|) p is exchanged with m. For reasons which will 



a) 2 -Aa/3 2 (1.10) 



become apparent immediately we have checked which solutions are stable with respect 
to the mean-field dynamics. 

In mean-field theory the transitions between the phases A and B respectively B and 
C are continuous and the transition between C and D is first-order. 

An exact calculation ]7| using the matrix product ansatz for the wave function 
on the line (3 = 1 has confirmed the existence of phases A and B. Limited Monte 
Carlo simulations |7j showed a strong similarity between the true and mean-field phase 
diagrams for a = 1. Since for the two-species model |2], [| [|, |9], §] the phase diagram 
obtained by mean- field is exact, it might be tempting to conclude that in the present 
case the mean-field phase diagram is exact. 

In Section 2 we show a detailed Monte Carlo analysis of the model with an 
unexpected result. The D phase exists only up to a line going through the diamonds 
in Fig.|l] where a first-order phase transition takes place. We will denote the domain 
under this line by D. The C phase does not exist, the B phase extends down to the 
D phase and will be denoted by B (see Fig.p^). In this picture the role of mean-field 
is quite perverse. In the D domain (broken phase) one has the two solutions described 
by equation ( |1.10p . Above the line of diamonds the system picks up the symmetric 
unstable solution given by equation (|1.8|) . This solution coincides with the one which 
defines phase B. 
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A useful tool in this analysis is played by the free energy functional which is defined 
as follows [|l(J: Take an order parameter w and let P(L,w) be the probability to have 
w for a lattice with L sites. Defining 

h(w) = -jhgP(L,w) (1.11) 

the free energy functional (FEF) is 

f(w) = lim f L (w) . (1.12) 

For equilibrium systems this is a convex function. As we are going to see this is not 
true for steady states which do not describe Gibbs ensembles. In the D domain if we 
take d as order parameter the FEF looks like a Ginzburg-Landau picture in the case of 



spontaneous breaking of a symmetry [fT9|| . On the first-order transition line separating 
the phases B and D the FEF is convex, looking like the FEF corresponding to a first- 
order phase transition in equilibrium phenomena. Above the diamond-line the FEF is 
a convex function with one minimum only. We have also looked at the spectrum of 
the hamiltonian (see Appendix B) on the separation line. It is gapless with a critical 



dynamical exponent |Tl| z = 2. This suggests, from the experience with the one-species 



model, the existence of shocks |T2], [13| . The problems of convergence coming from the 
definition (|1.12j ) are discussed in Appendix A. (Some of the ideas developed here are 



already in the paper of Bennett and Grinstein Jl4|). 

In Section 3 we study in detail (for a = 1) the nature of the first-order phase 
transition. A novel feature comes from the fact that we have not one but two order 
parameters (p and m or s and d). In the p-m-plane the minima of the FEF lie on 
a boomerang like figure. We show how this figure can be obtained by an interesting 
combination of shocks. The existence of shocks allows the prediction of different density 
profiles which are indeed observed in several Monte Carlo simulations. 

In Section 4 we consider the effect of a symmetry breaking term on the process (we 
take h bigger than zero in Eqs. (|1.2| - |1.3|) ). We show that the FEFs look like text book 
Ginzburg-Landau pictures including the existence of spinodal points in the D phase. As 



is well known jj5| such points do not exist in equilibrium phenomena. The mean-field 
calculations are presented elsewhere [|T6|. They are compared here with the Monte Carlo 
simulations. It looks like mean-field theory gives a good approximation to the data only 
for small values of (3 and large values of h. 

Moreover, we show that the spectrum of the hamiltonian defined in Appendix B 
shows massless exitations at the spinodal point with a dynamical critical exponent z = 1. 
This is a very surprising result for stochastic processes. This exponent is typical for 
continuous equilibrium phase transitions. In Section 5 we present our conclusions. 



2. The phase transition between the broken and unbroken phases 



In order to clarify the phase structure of the model, let us take a = 1. (This is a vertical 
line in Fig.|TJ) . If one wants to compare equilibrium with non-equilibrium phenomena, 
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it is useful to have in mind the two-dimensional n-state Potts model with n > 4 and 
interpret the parameter j3 as a temperature. Mean-field gives the following values of (3 
for the separations among the various phases: 

• f3 = 0.3289 (between phase D and phase C) 

• (3 = 0.3333 = 1/3 (between phase C and phase B) 

• (3 = 1 (between phase B and phase A) 

In order to check the mean-field predictions we have done Monte Carlo simulations to get 
fh(d) taking d (the difference between the averages of positive and negative particles) 
as order parameter. 

In Fig.^ we show our results for L = 400 and different values of (3 (all belonging to 
the D phase). In order to present the data in a clearer way, we have substracted the 
values of /l(0) from fi{d)- The data shown in Fig.[2] are interesting for two reasons: 
the shape of the FEF and the failure of mean-field predictions. First we notice that a 
dramatic change occurs around (3 = 0.274. As will be shown below the precise value 
is /3 cr it = 0.275. Below (3 CT \t the FEF is not a convex function, it has the behavior of 
a text-book Ginzburg-Landau picture, and shows two minima at the values given by 
mean-field (see Eq. ( |1.10| ) and the definition ( |1.5| )). 

At (3 cr \t the FEF behaves like the FEF of an equilibrium first-order phase transition. 
For (3 > (3 CV it the FEF has the expected behavior for a disordered phase. Notice that all 
the changes in the behavior of /l(^) occur at values of f3 below the value 0.3289 where, 



h(d) 



0.01 



0.005 



-0.005 




0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 A O.i 



Figure 2. f L (d) for =0.264, 0.268, 0.272, 0.274, 0.276, 0.280, 0.285, 0.290 and 
L = 400 sites. The values of /l(c?) are shifted by their value /l(0). 
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Figure 3. Free energy functional for (3 = 0.3: f L {d) for L =200, 400, 600 and 
1000 sites. The values of /l(cZ) for different L are shifted by /l(0). The little vertical 
stretches of the curves are due to the limited CPU time of the Monte Carlo runs. 



according to mean-field, the transition between the C and D phases is supposed to take 
place. 

In order to make sure that we are not dealing with finite-size effects we have 
computed fi(d) for 8 = 0.3 (still below 0.3289) for various number of sites up to 
L = 1000. The data are shown in Fig.^| (we have substracted the value /z,(0)). One 
can see that /l(^) is practically independent of L (the L dependence is in /l(0) - 
see Appendix A). This implies that one does not have finite-size effects, and that 
fi{d) — /l(0) is a convex function with one minimum as one expects to have if one 
is in the "disordered" phase. 

We now present more data which show what goes wrong with the mean-field 
prediction. In Fig.|] we show the positions of the minima of /l(g?) for various numbers 
of sites and various values of 8. One can observe that below (3 crit the minima converge 
to their mean-field values obtained from Eq.( |1.10|) and that a dramatic change takes 
place at 8 CT \ t . Above /3 cr i t the CP symmetry is not broken anymore (d — 0). 

In Fig.[5] we show the minima of fi{s) as a function of f3 for various L. The upper 
curve is given by the average value of s for the two phases (low/high density). This 



value is obtained from Eq.( |1.10| ): 

s = + - - [(1 + of - 4a/?] (2.1) 
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Figure 4. Minima of f L (d) for L = 100(o), 200(+), 400(d), 600(x ) sites. The curve 
( — • — ) gives the mean- field prediction obtained from Eq.( 1.10| ). 



The lower curve corresponds to the symmetric solution 

otf3 

s = 

a + /3 

As one can see from Fig.|5], below /3 cr i t the data converge to the values of s given by 
Eq.( |2.1|) and above /3 cr jt to the values of s given by Eq^^). This is very interesting 



unstable below (3 = 1/3: 
(2.2) 



because it means that fluctuations turn the unstable solution into a stable one and one 
obtains a first-order phase transition. The best estimate for /3 crit comes from the "fixed 
point" seen in Fig.|, one gets /5 crit = 0.275 ± 0.001 for a = 1. By the same method but 
not with the same patience we have computed the values of /3 cr it for other values of a. 
The values are given in Table 1. 



Table 1. The critical points for 


various a. 


a 


Pcrit 


a 


/?crit 


0.125 


0.11(2) 


1.00 


0.275(1) 


0.25 


0.186(2) 


1.25 


0.269(2) 


0.375 


0.234(2) 


1.50 


0.264(2) 


0.50 


0.259(2) 


1.75 


0.258(2) 


0.75 


0.274(2) 


2.00 


0.253(2) 
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Figure 5. Minima of ft(s) (s = (p + m)/2) as a function of /3 for L = 
100(o),200(+), 400(D), 600(x) sites. The upper curve corresponds to Eq.(fEl"l) (broken 
phase), the lower one to Eq.( |2.2| ) (unbroken phase). 



In the correct phase diagram which replaces Fig.l one has only three phases: the 
phase A like in Fig.|l], the phase B (the unbroken phase) which extends down to the 
diamond-line in FigJI] and finally the phase D under the diamond-line (see Fig.|TJ]). We 
have not checked if the seperation between the A and B phases is properly given by 
mean-field but we believe this to be the case because of the analytical calculation on 
the (3 = 1 line done in 0. 

A final test of our picture can be obtained looking at the spectrum of the hamiltonian 
(explicitly given in Appendix B). We consider the first three levels and take a = 1. Since 
the ground state has energy zero, the values of the energies of the next two levels, denoted 
by Ei and E 2 , coincide with the energy gaps. If our picture is correct Ei should vanish 
for large values of L for (3 < (3 CT \ t since in the large L limit one needs two states, the 
ground state and the first excited state, to describe the two vacua. E 2 should be finite. 
For (3 > (3 CT it both Ei and E 2 should be finite. Using the modified Arnoldi algorithm 
of Ref. [17] we have calculated, for various values of (3 and chain lengths (up to L = 11) 
the values of E x and E 2 . In order to get their large L limit we have used Bulirsch-Stoer 



approximants [ 18(1 - The results are shown in Fig.|6| which confirms the existence of a 



critical point around f3 = 0.275. 
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Figure 6. Estimates of the first (o) and second (+) excitations of the CP-symmetric 



hamiltonian (Eq. (7.1) with h = 0). The estimates were obtained from the spectra of 
the hamiltonian up to L = 11 by extrapolating the results in standard ways. Since all 
eigenvalues of H are positive their errors are estimated to be 0.0003. 



If one takes (3 = (3 cr - lt one can estimate the dynamical critical exponent z denned as 
follows: 

lim L Z E X {L) = const. > (2.3) 

L— *oo 

where E\(L) represents the energy of the first exited state for a chain of length L. We 
find z = 2.0 ± 0.1. This result is also interesting because from our experience with the 
one-species problem (see Appendix A) a value of z = 2 suggests the existence of shocks. 
In the next section we are going to find them. 

3. Description of the first-order phase transition 

We are now going to study in detail the point a = 1, (3 = /3 cr i t = 0.275 where as 
explained in the last section we have seen a first-order phase transition. The physics 
of the steady state can be understood only if one considers the two order parameters p 
and m. 

In Fig.[7| we show, in the p-m-plane, the contour lines of the FEF /l(p, m) near its 
minima. The figure suggests a large L limit in which we get the following picture: one 
line segment parallel to the p-axis (from U to V in Fig.0), one line segment parallel to 
the m-axis (from S to R), a smooth curve joining the points S, T and U. Let us give the 
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Figure 7. Contour lines of / L (p, m) for /3 = 0.275 and L = 1000 sites. The "height" 
difference between lines is constant. 



(p, m) coordinates of these five points: V=(0. 725, 0.149) where p and m are obtained from 
the mean-field equation ( |1.10| ), R= (0.149,0. 725), obtained from Eq. ( |1.10| ) exchanging p 
and m, T= (0.216,0.216), obtained from Eq. (|1.8|) . It is convenient to denote by pa and 
rriA the p and m coordinates of the point A. With this notation the coordinates of the 
points U and S are U=(l — py,my) and S=(pr, 1 — mn). The points V and R are to 
be expected to play a special role since they represent the broken phase for a value of j3 
slightly smaller than /3 cri t- Also the point T is to be expected since it corresponds to the 
symmetric phase (/3 slightly higher than /? cr i t ). The points S and U are unexpected since 
they do not correspond to any phase. Their coordinates have however a remarkable 
property: pu = 1 — py and similarly wig = 1 — ijir. These relations are typical for 
shocks (see Appendix A and Ref.[|12j). The following simple model captures 

the physics of the phase transition. First let us introduce the scaling variable z = k/L, 
where k is the site variable and L the length of the lattice. The model will apply in the 
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Figure 8. The shocks appearing at the first-order phase transition. The solid and 
broken lines give the profile of positive and negative particles, respectively. 



limit where k and L are large. We assume that the following processes take place: 

• With a probability Pi we have shocks of uncorrected positive particles as described 
by Fig.[8]a. The probability to find a positive particle at the left of the point g with 
a density pu and to the right of the point g with a density py is independent of z. 
One has fronts for < g < 1 with equal probability (this corresponds to the U-V 
line segment). The negative particles are also uncorrelated and have a constant 
density my. 

• Also with a probability Pi we have shocks of negative particles (see Fig.|8|b) which 
are just the CP-reflected shocks of the positive particles (this corresponds to the 




Figure 9. Density profile of the negative particles for L = 50(d), 100(o), 200(+) at 
the phase transition (3 = 0.275. The straight line is given by Eq.( |3.4D . 
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Figure 10. Density profile of the negative particles for L = 50(o), 100(+), 200(d) at 
j3 = 0.26, below the phase transition. The line gives the mean-field prediction, i.e. the 



average of the densities of Eq.(l.lC) 



R-S line segment). 

• With a probability P 2 we have configurations (Fig.||c) of constant density 
(independent of z) corresponding to all pairs of values p and m along the line 
segment S-U, all with the same probability. 

Obviously 

2P 1 + P 2 = 1. (3.1) 

The assumption that we have a line segment between S and U is an approximation 
since this line segment intersects the p = m line in the point of coordinates 
(ipv +Pr)/2, (m$ + mv)/2) = (0. 212, 0.212) which does not coincide with T but is very 
close to it. We did not try to find a better description of the S-T-U curve because this 
would imply many more Monte Carlo data than we were able to collect. 

The probabilities Pi and P 2 are determined from the condition that m) has 
the same value everywhere on the "boomerang" . We get 

Pi = P 2 = Pv ~ Pu r = 0.39 (3.2) 
2(pv - my) 

P 2 = = 0.22 . (3.3) 

pv — m v 

Let us now check the model, in this way we can present the data in an organized 
way. First, we can can compute f(d). One obtains a function which is flat between 
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— (pv — m v) an d (pv — wiy). This is just what one sees in Fig.§ for /? crit . Next we 
compute m(z) and get 

771(2) = Pim v + ^P 2 (^v + "is) + ^i(w R - z(m R - m s )) = 0.387 - 0.175 z . (3.4) 

Obviously p(z) = m(l — z). In Fig.^ we show the density of negative particles as a 
function of z for various numbers of sites. One notices that the straight line given by 
Eq.flSOD fits the data. 

In order to make sure that the linear density distribution observed in Fig.pl is not a 



finite-size effect, in Fig.|l(] we show for (3 = 0.26 (below the critical point) the function 
m(z). One observes that the data converge (slowly) from below to the mean-field value 
m = 0.44 obtained taking the average of the two equations ( |1.1U| ). 



In Fig.[TT| we show m(z) for (3 = 0.3 (above /3 cr ;t). The data converge from above to 
the value m = 0.23 which is obtained from Eq. (|1.8| ). In conclusion, the non-constant 
distributions seen just below and just above /3 cri t are cross-over phenomena. 

In order to "see" the shocks, we have looked at configurations where the density of 
positive particles p, averaged over the lattice, is fixed. This corresponds to a point on 
the U-V line segment in Fig.[7] and one expects to see a shock like in Fig.|8]a. Since we 
are also interested to see what kind of finite-size effects exist (see Appendix A), we have 
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Figure 12. Scaling of the shock: Density profiles of positive and negative particles 
for p = 0.5 fixed. Data for L =100, 200, 400 sites. The solid curve fitted to the shock 



is given by Eq.(3.7). Eq.(p.6[) gives the expression of ko- 
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Figure 13. Density profile p(z) and m(z) for d = 0.05 and L = 200. 
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chosen to present the densities as functions not of z but in terms of the variable 

V = (3.5) 

where ko is the position of the front which in the large L limit is given by 

k o= n L r (p-l + &rit)- (3.6) 

^Pcrit — J- 



The data for p = 0.5 are shown in Fig.|12[ Let us leave aside for a moment the existence 
of "tails" on the left hand side of the figure. The data scale, they are independent of L. 
A good fit to the data is given by 

V{y) = (\ - M er%/2) + I (3.7) 

Obviously, if we show the data as a function of z one gets the front from Fig.[8]a with 
g = 0.5. The value of m in Fig.[12| corresponds to the value of my in Fig.|8|a as expected. 
We now turn our attention to the "tails" . For the positive particles one can fit the data 
like 

p(k) = /xA^+/3 crit (3.8) 

where k is the position of the site starting with the left hand side of the chain and 
H = 0.08 is a constant. We have checked that for other shocks (one takes p different 
of 0.5 but still on the U-V segment of Fig.0) the value of /i stays unchanged. For 
negative particles the tails are exponential with a correlation length independent of L. 
This observation is intriguing because common sense would suggest that in this case 
one should have a finite correlation length also in the time direction, i.e. masses in the 
spectrum of the problem and not only the massless excitations discussed in Sec. 2. 

Another test of the model, this time on a part of it on which we certainly expect to 
be less precise, is to "cut" through the handle of the boomerang. We have taken d = 0.05 
in order to check if profiles like in Fig.[S]c are seen. The data are shown in Fig.|T^ together 
with the expected constant values (take the intersection of the line d = 0.05 with the 
line segment S-U in Fig.[7|). As one notices we find agreement although some mini-shocks 
are not excluded. After presenting the data we are left with a puzzle. We have expected 
the point T to be singled out in some way (it represents the symmetric phase), but it is 
not. 



4. The broken phase in the presence of a symmetry breaking external 
source 



As explained in Sec. 2 and shown in FigJTJ], the phase diagram of the CP-symmetric 



model consists of a power law phase A and a broken phase D separated by a first-order 
phase transition line from the unbroken phase B. In equilibrium systems, if one has a 
transition from an ordered to a disorderd phase, the broken phase corresponds to first- 
order phase transitions (and the FEF is a flat function between the two phases). As we 
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Figure 14. Phase diagram of the two-species model for h = obtained from the 
Monte Carlo simulations. The phase A is given by mean field. 



have seen this is not the case for steady states and for this reason we use the expression 
"first-order" only on the separation line where the FEF is flat and the expression "broken 
phase" for the domain where the FEF is not convex and has two minima. 

We consider now the effect of an explicit CP symmetry breaking in the model. We 
take h > in the boundary rates given by equations ( |1.2| ) and (|1.3|) . Using Monte Carlo 
simulations, we have computed the FEF /z,(<i, h) taking a = 1 and several values of (5 
and h] d is again the difference between the average densities of positive and negative 



particles. In Fig.|l5] we show fi(d,h) for (3 = 0.1 (in the D phase) and several values 
of h. We observe that below h crit = 0.2281 the FEF shows two minima and only one 
minimum for h > h cr \t. For h = h CT i t the first and second derivative of fL{d,h) with 



respect to d vanishes. This is the definition of a spinodal point [I5|. 

The behavior of fi{d,h) is the one expected from a Ginzburg-Landau ansatz [19 
Our data are of course not precise enough to give an estimate with four digits for /i crit , 
the value we give is obtained in mean-field. We will return soon to the problem of the 
validity of mean-field. 

In the Ginzburg-Landau ansatz f(d, h) has a simple h dependence: 

f(d,h)= f(d,h = 0) + vhd + g(h) (4.1) 



where v is a parameter. We have collected data for various values of (5 and h and find 
v = 1.4 getting a fit to the data within 10%. In Fig.IB we fix the value of h = 0.02 and 
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Figure 15. f L {d,h) for L = 30, a = 1.0, /3 = 0.1 and h = 0, 0.1, 0.2281, 0.3. For 
/i = the system is in the phase D. 



show the difference 

f L (d, h = 0) - f L (d,h = 0.02) (4.2) 



as a function of d. According to Eq.( [4.1| ) we expect straight lines with the slope —0.028. 
The data are compatible with this prediction. We have repeated the simulations for 
several values of h with the same result. 

We have computed the mean-field approximation of the model ( |l.l| - fL3|) for h 



different from zero [T^] , extending the calculation done in Ref . |7|] . Our motivation was 
to find out if, like in the h = case, one can understand in a simple way where and 
in which way mean- field is exact or wrong. In Table 2, for different values of (3 and 
h, we give the values of <i™ in which denotes the minima of the FEF as determined by 
Monte Carlo simulations as well as the values of c? m _f obtained in mean-field. Comparing 
the Monte-Carlo simulations with the mean-field results we reach the conclusion that 
probably mean-field is not exact anywhere. It is however a good approximation for small 
values of (3 and large values of h. 

We are now going to discuss some aspects of the dynamics of the model when one 
includes the symmetry breaking term. Our attention goes first to the spinodal point and 
we look at the spectrum of the hamiltonian (see Appendix B). We notice that the first 
excited state has an energy which vanishes in the large L limit. The dynamical critical 
exponent (see Eq.( [2.3|) ) determined using the values of the first excited state of chains up 
to 11 sites (3 = 0.1 and h CT i t = 0.22810 as given by mean-field is z = 1.00±0.01. We have 
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Figure 16. f L (d, h = 0) - f L (d, h = 0.02) for (3 =0.20, 0.26, 0.28, 0.30 and L = 100 
sites. 



repeated the same calculation for [3 = 0.05 (h crit = 0.28096) and found z = 1.000±0.001. 
As discussed previously, it is not clear if the values of h crit obtained by mean-field can 
be trusted to a precision of five digits. To our knowledge it is for the first time that the 
exponent z = 1 appears in stochastic processes. 

In Ref.|]20[| the toy-model of Godreche et al [^T| is examined (which corresponds to 
the (3 —>■ limit of the present model) and again z — 1 is found. More than that, in 
Ref. pOf it is shown that the spectrum has massless excitations which are universal and 
also massive excitations. (As opposed to the present model where we could study 
properly only one level, in the toy-model we can study many of them.) The toy- 
model has two parameters, one which corresponds to a in the original model and h. 
Normalizing the hamiltonian such that the sound velocity is the same, at the spinodal 
point the spectra corresponding to the massless regime do not depend on the former 
parameter. One can ask which implications has a z = 1 exponent. In Ref.pO| several 
properties of the time dependence of the order parameter are given. The simplest one 
being that the flip time from a configuration with the d value of the spinodal point 
to a configuration corresponding to the minimum of fi{d,h) (see Fig. |T5|) should be 
proportional to L. 

We would like to make another observation on the flip times related to barriers in 
the FEF. Let us take (3 = 0.1 and h = 0.1 in Fig. 15 and denote by Ai (A 2 ) the difference 
between the maximum of the FEF and the minimum of the FEF for negative (positive) 
values of d (A 2 = at the spinodal point). From the Monte Carlo data (L = 50 sites) 
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Table 2. The minima d™ n of /l(c£) (error ±0.01) and the mean-field predictions d m _f 
for a — 1, various /3 and /i. Values of d™ m predicted (within 10%) are printed bold 
face. Several values given for a set of parameters a, /3, h reflect different minima (d™ m ) 
or different stationary solutions of the mean-field equations (d m -f)- 
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0.6 





200 


0.00 
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200 


-0.08 


-0.0943 




0.35 





400 


0.00 









0.02 


100 


-0.05 


-0.468 






0.1 


100 


-0.45 


-0.512 




0.329 





400 


0.00 


-0.490/- 


hO.490/0 




0.02 


200 


-0.06 


-0.501 






0.1 


200 


-0.52 


-0.543 




0.3 





400 


0.00 


-0.537/- 


hO.537/0 




0.02 


400 


-0.15 


-0.546/+0.527/+0.098 




0.1 


400 


-0.58 


-0.584 




0.26 





200 


-0.53/+0.53 


-0.600/- 


1-0.600/0 




0.01 


100 


-O.53/+0.43 


-0.604/- 


h0.596/+0.019 




0.02 


100 


-0.54 


-0.607/- 


H0.592/+0.038 




0.07 


100 


-0.61 


-0.629 




0.1 





200 


-0.84/+0.84 


-0.849/- 


hO.849/0 




0.02 


100 


-0.84/+0.84 


-0.852/- 


H0.846/+0.00519 




0.05 


100 


-0.85/+0.83 


-0.856/+0.841/+0.013 




0.1 


50 


-0.86/+0.81 


-0.864/- 


H0.833/+0.026 



we get Ax = 0.178 ± 0.002 and A 2 = 0.041 ± 0.002. The two flip times from the stable 
(unstable) to the unstable (stable) phase are Ti ong (T short ). They have been measured 
for chains of different lengths L and their values are shown in Fig.[T7|. The straight lines 
in Fig.|l7| are given by the ansatz 

T long oc exp(Ai L) 

T short cx exp(A 2 L) . (4.3) 



The agreement between the ansatz ( |4.3| ) and the data shows a connection between the 
dynamics of the model (flip times) and the properties of the steady state (the FEF). 
Repeating the same procedure for a different value of j3 gives the same result. More 



about this subject can be found in Ref. [20 . 
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Figure 17. Logarithmic plot of the flip times for a — 1.0, (3 — 0.1 and h = 0.1. The 
dotted curves are given by 



5. Conclusions 

We have shown in the framework of the two-species model the following pattern of 
spontaneous CP symmetry breaking for steady states. If we start in the "disordered" 
phase (where the free energy functional as a function of the order parameter has only one 
minimum) and lower the "temperature" (the parameter (3) we have a first-order phase 
transition at /3 cr i t with a flat free energy functional (like in equilibrium problems). At /? cr it 
the density profiles and correlation functions in the scaling regime can be explained using 
a certain combination of shocks. Below the critical "temperature" /3 cr i t , in the broken 
phase, the free energy functional is not a convex function of the order parameter (unlike 
in equilibrium problems) but has two minima corresponding to two phases defined by 
their average densities of positive and negative particles. If the CP symmetry of the 
model is broken by a small perturbation (parameter h), one has only one minimum for 
the first-order phase transition (like in equilibrium) but still two minima in the broken 
phase. If we increase h one gets spinodal points. 

At the first-order phase transition we have seen massless excitations in the spectrum 
of the quantum chain hamiltonian associated with the dynamics of the process with a 
critical dynamical exponent z = 2. This does not exclude the presence of massive 
excitations which are suggested by the exponential tails seen in Fig.|D| and Fig{T^ in the 
density profiles. Massive excitations are not seen for the one-species model described in 
Appendix A [f£|. 
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At the spinodal points we have also seen massless excitations with z = 1 as well as 
massive excitations. The latter are suggested by exponential tails observed in the density 
profiles not shown in this paper. These observations are confirmed by the investigation 
of the low (3 limit of the model |2(J . At this point it is not clear if the exponent z = 1 has 
profound implications like in equilibrium problems where it implies conformal invariance 
or is just one dynamical critical exponent among others. We have looked to conformal 
towers in the spectrum of the hamiltonian but the data of our short quantum chains 
did not give any clean results. This is probably due to the fact that the critical point 
is not known exactly (we used the values obtained from mean-field) and that in general 
the convergence to the thermodynamic limits is slower in stochastic quantities than in 
equilibrium ones. This is an empirical observation. We think that a further investigation 
of the model is necessary by taking the second rate in ( |1 . 1[ ) not equal to the other two. 
Some preliminary results are described in Ref.fTj. This study can be done not only by 
Monte Carlo simulations but also using the analytical methods developed recently in 
Ref.lEl. 
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6. The asymmetric exclusion process 

The physics of this model is well understood [0, here we are going to discuss some 
aspects of it which are necessary for the understanding of the finite-size effects of the 
two-species model. We consider a chain of length L with two types of particles (+ and 
— ). In an infinitesimal time step dt the following processes take place: 

(+)* (-)fe+i -> (-)fc (+)fc+i with probability dt 

(— )i — > (+)i with probability a dt (6.1) 

— > (— )l with probability (3 dt. 



Defining o = a + (3 and h = (3 — a one notices that for h = the model is CP- 
invariant. We denote by c the average density of positive particles. The model presents 
three phases: 

• the power law phase (a > 1/2, (3 > 1/2) where c = 1/2 

• the low density phase (a < 1/2, (3 > a) where c = a 

• the high density phase (a > /5, (3 < 1/2 ) where c = 1 — (3 

A coexistence line (a = (3 < 1/2) separates the low density from the high density phase. 
The concentration has a discontinuity along this line from c = /3 to c=l — (3. 

Using Monte Carlo simulations we have computed the FEF /l(c, h) keeping a fixed 
and using several values of h. The results are shown in Fig.lTB]. We notice that for h = 0, 
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Figure 18. Asymmetric diffusion model: /l(c, h) for s = 0.25 and h = 0, 0.004, 0.2, 
each for L = 800 and 1000 sites. We subtracted /i(c m in, h) from each curve. 



the FEF is a flat function between c = (3 and c = 1 — (3. A small positive h (h = 0.004) 
gives a FEF with only one minimum at c = (3 just like in equilibrium first-order phase 
transitions. A larger value of h (h = 0.1) moves the minimum away. We have studied 
the convergence of the function fi(c, h) to f(c, h) taking c m i n the value of c where the 
FEF has a minimum. For all the examples we have studied we have found that 

/i(c m in, h) = f(c min , h) + a/L (6.2) 

Where the constant a depends on o and h. For example if we take a = 0.5 and h = 0.2 



like in Fig.|18|we find a = 4.2±0.2. We have also found that if we start on the coexistence 
line with a certain value of a and switch-on the CP symmetry breaking term h, we have 
the following simple expression for the FEF: 

f(c,h)=f(c,h = 0) + uch + g(h) (6.3) 

where v depends on a. For the example a = 0.5 studied above we find v = 4.5. We 
have observed that if we take L = 200 the FEF has already a shape independent of L, 
the single L dependence being in /^(cmin, h). 

We would like to mention that much of the physics which will be described below 



is contained in a simplified model studied by Schiitz [24]. In this model the bulk 
dynamics is deterministic and only the boundaries are stochastic. The critical dynamical 
exponent on the coexistence line of the present model is z = 2 like in Ref . |24[] . This 



corresponds to the picture in which the front of a shock performs a symmetric random 
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c(y) 




Figure 19. Scaling of the shock: Density profiles for a = j3 = 0.25 (with c = 0.5 
fixed) and for L =100, 200. The dotted curve is 0.25tanh(y) + 0.5. The wriggled 
curves represent the Monte Carlo data which coincide with the solid curve within 1%. 



walk. 

On the coexistence line the density of positive particles is high at the right hand side 
of the chain (c = 1 — 0) and low at the left hand side of the chain (c = /?,/?< 1/2). The 
two regimes are joined by a shock front. One measures such a shock by constraining 
the average density of positive particles to a given value between (3 and 1 — (3. If c(k) is 
the concentration at lattice site k, it is useful to define the variable 



y 



k — k 

~7T 



(6.4) 



where 



kn 



2(3-1 



1 + P) 



(6.5) 



is the position of the shock front. We have found empirically that the Monte Carlo data 
are fitted (within 1%) by the function 



c(y) = (i - /3) tanh(y) + i . 



(6.6) 



In Fig.fBl we show the Monte Carlo data for c = 0.5 and a — (3 — 0.25. The existence 



of a front with width of order L 1 / 2 is to be expected from the model studied by Schiitz 
2M , from the exclusion model with a blockage [EH B6J and from the asymmetric model 
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with an impurity [27]. Our purpose in studying the finite-size behavior of the shocks 
was to find out how many lattice sites are necessary to be in the scaling limit in order to 
know for the three-states model which lattice size are necessary to see shocks. It turns 
out that L = 100 are enough. 



7. The Hamilton operator 



The time evolution operator H of the master equation (|1.4j) for the processes fll.lQ - 
( |1.3|) can be written as a sum of operators acting non-trivially on adjacent sites only 
(hopping) and of two operators acting on the first and last site (input and output) M: 



L-l 



H = Hi + H kjk+ i + H L 

k=l 



The matrices are 



H, 



k,k+l 



h 



Ik-1 



H L = h 



a 
-a 




lL-1 



(7.1) 
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where the subscripts denotes the sites the matrices act on and I is the 3x3 unit matrix. 



The order of the basis at each site is chosen as 
hermitian. 



>, |+ >, |— >. Note that H is not 
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